﻿#include <iostream>
#include <fstream>
#include <math.h>
#include <exception>
#include <string>
#include <vector>
#include <sstream>
#include "func.h"
#include "parm.h"


using namespace std;

int main(int argc, char** argv)
{
	//读取a0,th2,rin,rout,积分距离r2,团块环向宽度dphcl,团块个数t_num,输入文件名filename_in,输出文件名filename_out
	//输出 理论上qpo的中心频率，N*sum（x-x_i）^2,N_ph即直流分量 
	//此时的N_ph还应该加上其余辐射情况 除sum才等于 真正的rms
	double a0 = stod(argv[1]);
	double th2 = stod(argv[2]);
	double rin = stod(argv[3]);
	double rout = stod(argv[4]);
	double r2 = stod(argv[5]);
	double dphcl = stod(argv[6]);
	int t_num = stoi(argv[7]);

	ifstream input(argv[8]);

	double z1 = 1 + pow(1 - a0 * a0, 1.0 / 3.0) * (pow(1 + a0, 1.0 / 3.0) + pow(1 - a0, 1.0 / 3.0));
	double z2 = sqrt(3 * a0 * a0 + z1 * z1);
	double r_min = 3 + z2 - sqrt((3 - z1) * (3 + z1 + 2 * z2)) + 0.001;



	double rcl0 = 0.5 * (rin + rout);

	//计算相关变量

	double omg = omega(a0, rcl0);
	double period = 2 * Pi / omg;
	double dt = period / time_resolution;
	double rclp = pow(rcl0, p0);
	dphcl = dphcl * Pi / 180;

	double flux[time_resolution] = {};

	try
	{
		while (true)
		{
			double data[data_len] = {};

			for (int i = 0; i < data_len; ++i) {
				if (!(input >> data[i])) throw runtime_error("read over");
			}

			const double phi1 = phiem(data[xt], data[yt], data[sr], data[stheta]);
			//原团块对应的时间
			double time_origin_end = phi1 / omg;


			//像对应的时间
			double time_image_end = time_origin_end + data[str] + data[stth] - r2;
			double time_image_begin = time_image_end - dphcl / omg;
			int time_image_begin_int = normalization_time(time_image_begin, period, time_resolution);
			int time_image_end_int = normalization_time(time_image_end, period, time_resolution);
			//计算能段和强度
			//int ji = round(elin * data[gt] / 0.1);
			double sg = data[ds] * pow(data[gt], 3) / rclp;

			//一周期均匀分布的团块数t_num，因此两个团块相隔的单位时间数为，1000/t_num

			for (int t_count = 0; t_count < t_num; t_count++)
			{
				int del_t = t_count * time_resolution / t_num;
				for (int t_image = time_image_begin_int + del_t; t_image < time_image_end_int + del_t; t_image++)
				{
					int t_in_period = t_image % time_resolution;
					if (t_in_period < 0)
					{
						t_in_period += time_resolution;
					}
					flux[t_in_period] += sg;
				}
			}
			//将符合时间要求(在第一个周期)的光子 加权加入结果中
		}
	}
	catch (runtime_error err)
	{
		input.close();
	}
	double sum_square = 0, N_ph = 0;
	fractional_rms(flux, time_resolution, &sum_square, &N_ph);
	ofstream output(argv[9]); //最终输出的二维函数
	output << period << endl;
	output << sqrt(sum_square * time_resolution) << endl;//该值与N_ph的比值为farctional rms
	output << N_ph << endl;//后续应该再与外盘的流量相加才是真正的N_ph
	output.close();

	return 0;
}

